#include "mex.h"
#include "cc.h"

void run_linkage(float *output, float *input, unsigned int size_y);

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	double real_size_y;
	double *inMatrix;
	unsigned int size_x, size_y;
	double *outMatrix;
	float *in, *out;

	if (nrhs != 1) {
		mexErrMsgTxt("One input required.");
	}

	if (nlhs > 1) {
		mexErrMsgTxt("Cannot have more than one output.");
	}

	/*
	if (mxGetM(prhs[0]) == 1 && mxGetN(prhs[0]) == 1) {
		mexErrMsgIdAndTxt("CudaCluster:linkage:notRowVector",
				"Input must be a matrix.");
	}
	*/

	inMatrix = mxGetPr(prhs[0]);
	size_x = mxGetN(prhs[0]);
	size_y = mxGetM(prhs[0]);

	// solve for y:  0 = y^2 - y - 2*size_x
	real_size_y = (1 + sqrt(1 + 8*(int)size_x))/2;

	out = mxMalloc(sizeof(float)*(real_size_y-1)*3);
	if (!out)
		mexErrMsgTxt("Could not alloc out.");
	in = mxMalloc(sizeof(float)*size_x*size_y);
	if (!in) {
		mxFree(out);
		mexErrMsgTxt("Could not alloc in.");
	}
	plhs[0] = mxCreateDoubleMatrix(real_size_y-1, 3, mxREAL);
	outMatrix = mxGetPr(plhs[0]);
	if (!outMatrix) {
		mxFree(out);
		mxFree(in);
		mexErrMsgTxt("Could not alloc outMatrix.");
	}

	double2float(in, inMatrix, size_x, size_y, 0);

	run_linkage(out, in, real_size_y);

	float2double(outMatrix, out, 3, real_size_y-1, 1);

	mxFree(in);
	mxFree(out);
}


